/*
 * solver.c
 *
 *  Created on: 30/12/2011
 *      Author: marcio
 */

#include "solver.h"

extern long k;
extern long* va;

double** x;
double* mu;
double* subpb;
int noI;//melhoria
double bestLB;//o melhor limite inferior encontrado
double bestUB;//o melhor limite superior encontrado
double lb;//o limite inferior corrente
double ub;//o limite superior corrente
double pi;//passo do lagrangeano

int compar(const void* a, const void* b){
	int v = *(int*)a;
	int u = *(int*)b;
	return mu[v] - mu[u];
}


double calculateUB(){
	setAdjMatrix(getSharedGraph());
	int V[nverts()];
	int cor[nverts()];
	int i,j;
	int color = 1;
	int adj[nverts()];
	int cand = -1;
	for(i = 0; i < nverts(); i++){
		V[i] = i;
		cor[i] = 0;
		adj[i] = 0;
	}
	qsort(V, nverts(), sizeof(int), compar);
	cor[0] = 0;
	for(i = 1; i< nverts(); i++){
		for(j = 0; j < i; j++){
			if(!hasEdge(V[i], V[j])){
				adj[cor[j]] = 1;
			}
		}
		for(j = color - 1; j >= 0; j--){
			if(adj[j] == 0)cand = j;
		}
		if(cand == -1){
			cand = color;
			color++;
		}
		cor[i] = cand;
		cand = -1;
		for(j = 0; j < nverts(); j++){
			adj[j] = 0;
		}
	}
	if(bestUB > color){
		bestUB = color;
	}
	return color;
}

int noImproviment(){
	noI++;
	if(noI > 30){
		pi = pi/2.0;
		noI = 0;
	}
	return 0;
}

double calculateLB(){
	int i;
	double value  = 0;
	for(i = 0; i < nverts(); i++){
		//printf("value = %f \n", value);
		//printf("mu[i] = %f \n", mu[i]);
		//printf("subpb = %f \n", subpb[i]);
		value = value + mu[i] + subpb[i];
	}
	if(value > bestLB) {
		printf("Atualizei LB para %f\n", value);
		bestLB = value;
	}
	else noImproviment();
	//printf("lb = %f e bestLB = %f \n", value, bestLB);
	//getchar();
	return value;
}

int updateMuLambda(){
	setAdjMatrix(getSharedGraph());
	double g[nverts()];
	double theta = 0.0;
	double theta_div = 0.0;
	int i,j;
	for(i = 0; i < nverts(); i++){
		g[va[i]] = 1;
		for(j = 0; j <= i; j++){
			g[va[i]] -= x[va[j]][va[i]];
		}
		//printf("g[%lu] = %f \n", va[i], g[va[i]]);
		theta_div += g[va[i]]*g[va[i]];
	}
	//printf("theta_div= %f \n", theta_div);
	if(theta_div >=  0.001)theta = pi*(ub - lb)/theta_div;
	else return 1;
	//calcula subgradientes
	for(i = 0; i < nverts(); i++){
		mu[va[i]] =+ theta*g[va[i]];
		if(mu[va[i]] < 0.0000){
			mu[va[i]] = 0.0;
		}
		if(mu[va[i]] > 1.0){
			mu[va[i]] = 1.0;
		}
		//printf("mu[%lu] = %f \n", va[i], mu[va[i]]);
	}
	return 0;
}


/*cria a formulação*/
int createFormulation(){
	setAdjMatrix(getSharedGraph());
	int i, j, ckeck,m;
	pi = 2.0;
	bestLB = -nverts();
	bestUB = nverts();
	mu = malloc(sizeof(double)*nverts());
	subpb = malloc(sizeof(double)*nverts());
	x = malloc(sizeof(double*)*nverts());
	for(i = 0; i < nverts(); i++){
		mu[i] = 1.0;
		x[i] = malloc(sizeof(double)*nverts());
	}
	init();
	for(i = 0; i < nverts(); i++){
		for(j = i+1; j < nverts(); j++){
			if(hasEdge(va[i], va[j])){
				ckeck = 0;
				for(m = i+1; m < nverts(); m++){
					if(!hasEdge(va[j], va[m]) && hasEdge(va[i], va[m]) && m!= j){
						ckeck = 1;
						if(m > j)addEdgeConstrain(va[i],va[j],va[m]);
					}
				}
				if(ckeck == 0){
					addIndConstrain(va[i],va[j]);
				}
			}
		}
	}
	return 0;
}

int solveSubProblem(){
	int i;
	for(i = 0; i < nverts(); i++){
		subpb[i] = solveSingleLP(i);
		//printf("subp = %f \n", subpb[i]);
		//getchar();
	}
	return 0;
}

/*resolve a formulação*/
double solve(){
	while(pi > 0.005){
		changeObjCoef(mu);
		solveSubProblem();
		getX(x);
		lb = calculateLB();
		ub = calculateUB();
		//printf("lb = %f e ub = %f \n", lb, ub);
		if(updateMuLambda()!= 0) return bestLB;
	}
	return bestLB;
}
